!nan***********************************************************************
      subroutine nanChecker(jd,hr,val,idx)

!     tested with Intel- and g77-compilers
!--------------------------------------------------------------------------
      implicit none

      character*3 nancheck
      double precision val,hr
      integer jd,idx

      write(nancheck,'(f3.0)') val
      if (nancheck(1:1).eq.'?'.or.nancheck(1:1).eq.'n'.or. &
          nancheck(1:1).eq.'N'.or.nancheck(1:2).eq.' N') then
         if (idx.eq.1) val=1.9
         if (idx.eq.2) val=-50.
         print*,'nan in ',jd,hr,idx
      endif

      return
      end




!gfx++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine report(xpsi)
!
!     Calculates summaries and prints simulation results
!-----------------------------------------------------------------------
      implicit none

      common/res/transp,co2_fx,wdnlay,cantmp,cantotgl
      common/pcontr/alat,date,live
      common/plight/lai,sai,angst
      common/vligt1/aldif,ttsi,sarea
      common/senerg/par(2),tpairc,swtot(2),templf(2),therm,tran(2)
      common/vphoto/ci(2),ga,pn(2),gs(2),gl(2),wd(2)

      double precision tralay,psnlay,area(2),sunfr,shdfr,cantotgl
      double precision transp,co2_fx,cangl,cantmp,avggl,avgtmp,xpsi
      double precision alat,date,live,lai,sai,angst,aldif,ttsi,sarea, &
          wdnlay,par,tpairc,swtot,templf,therm,tran,ci,ga,pn,gs,gl,wd
      integer j


!     Calculate sunlit and shaded leaf area
!     LAI        Leave area index for layer L
      area(1)=sarea
      area(2)=lai-sarea

!     TRAN(J)       Transpiration on a leaf area basis   [mmol m-2 s-1]
!     AREA(J)       Area of sun, shade, or all leaves    [m2]
!     TRALAY        Transpiration on a ground area basis [moles m-2 s-1]
      tralay=0.0
      psnlay=0.0
      wdnlay=0.0
      cantmp=0.
      cantotgl=0.

!     h2o-coupling to soil layer
!      live=0.0142*xpsi+1.7149
!      live=0.0142*xpsi+1.6              ! ok for all
      live=0.0142*xpsi+1.7149
      if (live.lt.0.) live=0.
      if (live.gt.1.) live=1.

!     Consider both sunlit (J=1) and shaded (J=2) leaves
      do j=1,2
!        Transpiration of leaves in layer
!        (moles of H2O per square meter ground area per second)
!        Sum up layer totals for PS and TR
         tralay=tralay+tran(j)*area(j)*live

!        Calculate ps in layer on ground area basis
!        Units: umol * m-2 meter horizontal area * s-1
         psnlay=psnlay+pn(j)*area(j)*live

         wdnlay=wdnlay+wd(j)*area(j)*live
      enddo


!---- CANGL Average canopy conductance [mmol m-2 s-1] -------
!     Fractions of leaf area that's sunlit, shaded.
      if (lai.gt.0.) then
         sunfr=area(1)/lai
         shdfr=1.-sunfr
      else
         sunfr=0.
         shdfr=0.
      endif

!     Average leaf conductance, not calculated from transpiration,
!     per leaf area
      avggl=gl(1)*sunfr+gl(2)*shdfr

!     total canopy conductance [mmol m-2 s-1]
      if (lai.gt.0.) then
         cantotgl=avggl*lai
      else
         cantotgl=0.
      endif

!     Average Temperature
      avgtmp=templf(1)*sunfr+templf(2)*shdfr
      if (lai.gt.0.) then
         cantmp=avgtmp
      else
         cantmp=0.
      endif

!     total Transpiration
!     3.6 converts millimoles/s to moles/h
!     1.8e-2: mol/m2*h --> l/m2h
!     1.8e-5: mol/m2*h --> m/h=m3/m2
!     0.3 mm entspricht 206 W/m2
      transp=tralay*3.6*1.8e-2

!     CO2-Flux
!     umol/m2s      +wdnlay=GPP
!     sresp im umol m-2 s-1
!cc      ssresp=sresp*3.6
      co2_fx=psnlay+wdnlay


      return
      end


!spat===================================================================
      subroutine solar(longitude,rad_lat,julian,hour,slope,aspect, &
          cosThetaNormal,cosTheta,omega)
!
!-----------------------------------------------------------------------

      USE constants

      implicit none

      integer julian,hour
      double precision slope,aspect
      double precision rad_lat,longitude,declination,asp,slp
      double precision cosTheta,omega,Lo
      double precision cosThetaNormal       ! ,solan


!     W m-2
!     SolarConstant=1360.

      cosTheta=0.0
      asp=aspect
      slp=slope

!       N   E   S  W    N      (slope aspect)
!     -180 -90  0  90  180     (slope aspect)
!       0   6  12  18   24 h   time
      asp=asp-180.

!---------------------only valid for MEZ-----------------------
!     solar time in the time zone MEZ (15 degree longitude)
!     correct for deviations from the 15 degree meridian
!      omega=hour+(1.-longitude/360.*24.)
!
!     scale 0 -> 24 hours to  -180 -> 180
!     9:00 a.m. == -45 degrees
!      omega=omega*360./24.-180.
!--------------------------------------------------------------

!c     LINACRE,E.: Climate, Data and Resources
!c     hour=Local Standard Time (Tls)
!c     Tt=True Solar Time     Tu=Greenwich Time
!c     Tlm=Local Mean Time (Tt-E)
!c     E=Value of the Equation of Time (Smithsonian Meteorological Tables)
!c     L=Longitude of place
!c     Lo=Longitude of the standard meridian for the place,
!c        defining the local time zone
!c     Tt=Tu-L/15+E/60        Tls=Tu-Lo/15     omega=15*(Tt-12)
      Lo=int(longitude)/15*15
!c     omega=hour angle in degrees. Equation of Time (E) ignored (max 20min)
      omega=15.*((Lo/15.+hour-longitude/15.)-12)
      omega=omega*Pi/180.

!c     LINACRE: declination=23.45*sin(0.986*Nm*Pi/180.) in degrees
!c     Nm=number of days since 21 March    0.986=360/365
      declination=23.45*sin(Pi/180.*360.*(284.+julian)/365.)
      declination=declination*Pi/180.

!c     if slope is in degrees, uncomment the following line
      slp=slp*Pi/180.

      asp=asp*Pi/180.
!cgy      rad_lat=latitude*Pi/180.

!c      print*,'latitude: ',latitude
!c      print*,'declination: ',declination
!c      print*,'slope: ',slp
!c      print*,'aspect: ',asp
!c      print*,'omega: ',omega

      cosThetaNormal=sin(declination)*sin(rad_lat)+cos(declination)*  &
          cos(rad_lat)*cos(omega)

!c      print*,'Cosine of zenith angle on flat surface is ',cosThetaNormal
!c      print*,'In degrees:  ',acos(cosThetaNormal)*180/Pi

      cosTheta=(sin(rad_lat)*cos(slp)-cos(rad_lat)*sin(slp)* &
          cos(asp))*sin(declination)+(cos(rad_lat)*cos(slp)+ &
          sin(rad_lat)*sin(slp)*cos(asp))*cos(declination)* &
          cos(omega)+sin(slp)*sin(asp)*cos(declination)* &
          sin(omega)

!c      print*
!c      print*,'Cosine of zenith angle is ',cosTheta
!c      print*,'In degrees:  ',acos(cosTheta)*180/Pi
!c      print*,'Direct Irradiaton on slope: ',SolarConstant*cosTheta
!c      solan=Pi/2.-acos(cosTheta)
!c      print*,'Solar angle: ',solan
!c      print*,'In degrees: ',solan*180./Pi,slp*180./Pi,asp*180./Pi
!c      if (solan.gt.0.) read(*,*)

      return
      end







!soil--------------------------------------------------------------------
      subroutine readSoilParas(smatrix,rmatrix)

!     read soil inputfiles into matrix
!gy   besser alles in eine Schleife oder alles in eine Datei !
!c------------------------------------------------------------------------
     IMPLICIT NONE

      integer nos,nsp,nrp,nsl

      parameter (nos=1)         ! number of soils
      parameter (nsp=16)        ! number of soilparameters
      parameter (nrp=10)        ! number of soilrespparameters
      parameter (nsl=4)         ! number of soillayers + 1 !!!

      character*200 line
      integer status,li,lli
      double precision smatrix(nsp,nsl,nos),rmatrix(nrp,nsl,nos)

!c     Bodenparameter lesen
      open(13,file='eag_mult.bod',status='old')
      do li=1,nos
!c        Parameter fuer Schichten
         do lli=1,nsl
            call getline(13,line,status)
            if (status.eq.1) goto 20

!c           thick xtheta rootdens thetaS kSat vgN vgAlpha thetaR bd som xtemp lambA lambB lambC lambD lambE
!c             1      2      3        4     5   6     7       8   9   10   11    12    13    14    15    16
            read(line,*) smatrix(1,lli,li),smatrix(2,lli,li), &
                smatrix(3,lli,li),smatrix(4,lli,li), &
                smatrix(5,lli,li),smatrix(6,lli,li), &
                smatrix(7,lli,li),smatrix(8,lli,li), &
                smatrix(9,lli,li),smatrix(10,lli,li), &
                smatrix(11,lli,li),smatrix(12,lli,li), &
                smatrix(13,lli,li),smatrix(14,lli,li), &
                smatrix(15,lli,li),smatrix(16,lli,li)

            smatrix(5,lli,li)=smatrix(5,lli,li)*36000.      ! kSat is mm/h

            if (smatrix(1,lli,li).eq.0.) goto 10
         enddo
 10      continue
      enddo
 20   close(13)


!c     Bodenatmungsparameter lesen
      open(13,file='boden.atm',status='old')
      do li=1,nos
!c        Parameter fuer Schichten
         do lli=1,nsl
            call getline(13,line,status)
            if (status.eq.1) goto 40

!c           Corg k1 k2 fractk1 respTP1 respTP2 respTP3 respFP1 respFP2 respFP3
!c            1   2  3     4       5       6       7       8       9      10
            read(line,*) rmatrix(1,lli,li),rmatrix(2,lli,li), &
                rmatrix(3,lli,li),rmatrix(4,lli,li), &
                rmatrix(5,lli,li),rmatrix(6,lli,li), &
                rmatrix(7,lli,li),rmatrix(8,lli,li), &
                rmatrix(9,lli,li),rmatrix(10,lli,li)

            if (rmatrix(1,lli,li).lt.-9990.) goto 30
         enddo
 30      continue
      enddo
 40   close(13)

      return
      end





!csoil--------------------------------------------------------------------
      subroutine initStartValues(smatrix,rmatrix,bodart)

!c     give some common variables some startvalues
!c------------------------------------------------------------------------

      USE soil1pix_dble

      implicit none

      integer nsp,nrp

      parameter (nsp=16)        ! number of soilparameters
      parameter (nrp=10)        ! number of soilrespparameters


      real theta0,sst0
      double precision maxPsi,maxPsiForDt,thetaPsi,psiTheta
      double precision smatrix(nsp,nsl,nos),rmatrix(nrp,nsl,nos)
      integer i,k,bodart

      theta0=0.3
      sst0=3.

      do k=1,nos
         call initSoilCommon(smatrix,rmatrix,k)

!c        Berchnung des Wassergehalts bis zu dem der Boden maximal
!c        aufgefuellt werden darf, damit Zeitschritte nicht zu kurz werden
!c        userMaxPsi=-(10.0**2.5)
         maxPsi=maxPsiForDt(1./60.) !=1min

!c        auch der Benutzer kann ein maximales Psi vorgeben
!c        (z.B. Feldkapazitaet: Bucket!)
         maxPsi=min(maxPsi,-50.)         ! usermaxpsi=-50

         do i=1,noflayers
            thetaMax(i,k)=thetaPsi(maxPsi,i)
         enddo

         do i=1,noflayers
            if (k.eq.bodart) then
               theta(i)=min(theta0,thetaMax(i,k))
               temp(i)=sst0
               psi(i)=psiTheta(theta(i),i)
            endif
         enddo
      enddo

      return
      end






!soil--------------------------------------------------------------------
      subroutine initSoilCommon(smatrix,rmatrix,styp)

!c     init soil common-blocks
!c------------------------------------------------------------------------

      USE soil1pix_dble

      implicit none

      integer nsp,nrp

      parameter (nsp=16)        ! number of soilparameters
      parameter (nrp=10)        ! number of soilrespparameters


      integer styp,lli
      double precision xtheta(nsl),xtemp(nsl),zzz
      double precision smatrix(nsp,nsl,nos),rmatrix(nrp,nsl,nos)

!     Parameter fuer Schichten
      do lli=1,nsl
         Corg(lli)=rmatrix(1,lli,styp)
         k1(lli)=rmatrix(2,lli,styp)
         k2(lli)=rmatrix(3,lli,styp)
         fractk1(lli)=rmatrix(4,lli,styp)
         respTP1(lli)=rmatrix(5,lli,styp)
         respTP2(lli)=rmatrix(6,lli,styp)
         respTP3(lli)=rmatrix(7,lli,styp)
         respFP1(lli)=rmatrix(8,lli,styp)
         respFP2(lli)=rmatrix(9,lli,styp)
         respFP3(lli)=rmatrix(10,lli,styp)

         thick(lli)=smatrix(1,lli,styp)
         xtheta(lli)=smatrix(2,lli,styp)
         rootdens(lli)=smatrix(3,lli,styp)
         thetaS(lli)=smatrix(4,lli,styp)
         kSat(lli)=smatrix(5,lli,styp)
         vgN(lli)=smatrix(6,lli,styp)
         vgAlpha(lli)=smatrix(7,lli,styp)
         thetaR(lli)=smatrix(8,lli,styp)
         bd(lli)=smatrix(9,lli,styp)
         som(lli)=smatrix(10,lli,styp)
         xtemp(lli)=smatrix(11,lli,styp)
         lambA(lli)=smatrix(12,lli,styp)
         lambB(lli)=smatrix(13,lli,styp)
         lambC(lli)=smatrix(14,lli,styp)
         lambD(lli)=smatrix(15,lli,styp)
         lambE(lli)=smatrix(16,lli,styp)

         if (vgN(lli).gt.0) vgM(lli)=1.0-1.0/vgN(lli)

         if (thick(lli).eq.0.) then
            nOfLayers=lli-1
!           drainage (0) or water table (1)
            lowBouZero=int(xtheta(lli))
            goto 10
         endif
      enddo
 10   continue

      if (nOfLayers.lt.2) write(*,*)  &
          'Kopplung von gfac an cpsi(2) nicht moeglich!'

!c     Uebertragung auf neue Bodendatenstruktur
      zzz=0.0
      do lli=1,nOfLayers
         bb(lli)=-(3.0*(vgN(lli)-1.0)+2.0)
         zzz=zzz+thick(lli)
         zz(lli)=zzz
      enddo

      return
      end




!crop***********************************************************************
      subroutine readCropFiles(cmatrix)
!
!     read crop specific data from parameter file
!--------------------------------------------------------------------------
      implicit none

      integer cspec,ccols,crows

      parameter (cspec=1)
      parameter (ccols=6)
      parameter (crows=10)


!      integer i  ! commented by Bumsuk
      character*15 filename
      double precision cmatrix(ccols,crows,cspec)

      filename='rice_odae.crp'

      open(unit=95,file=filename,status='old')
      call readCropParas(95,cmatrix)
      close (unit=95)

      return
      end







!crop***********************************************************************
      subroutine readCropParas(fp,cmatrix)
!
!     read crop specific data from parameter file
!--------------------------------------------------------------------------

      implicit none

      integer cspec,ccols,crows

      parameter (cspec=1)       ! max. number of crop species
      parameter (ccols=6)       ! number of columns in the cropparameterfile
      parameter (crows=10)      ! number of rows in the cropparameterfile


      character*200 line
      integer fp,status,i
      double precision cmatrix(ccols,crows,cspec)

      status=0

      call getline(fp,line,status)
      if (status.eq.1) goto 20
      read(line,*) (cmatrix(i,1,1),i=1,3)

      call getline(fp,line,status)
      if (status.eq.1) goto 20
      read(line,*) (cmatrix(i,2,1),i=1,5)

      call getline(fp,line,status)
      if (status.eq.1) goto 20
      read(line,*) (cmatrix(i,3,1),i=1,4)

      call getline(fp,line,status)
      if (status.eq.1) goto 20
      read(line,*) (cmatrix(i,4,1),i=1,5)

      call getline(fp,line,status)
      if (status.eq.1) goto 20
      read(line,*) (cmatrix(i,5,1),i=1,5)

      call getline(fp,line,status)
      if (status.eq.1) goto 20
      read(line,*) (cmatrix(i,6,1),i=1,5)

      call getline(fp,line,status)
      if (status.eq.1) goto 20
      read(line,*) (cmatrix(i,7,1),i=1,5)

      call getline(fp,line,status)
      if (status.eq.1) goto 20
      read(line,*) (cmatrix(i,8,1),i=1,6)

      if (int(cmatrix(1,1,1)+0.5).eq.1) then  ! grassflag
         call getline(fp,line,status)
         if (status.eq.1) goto 20
         read(line,*) (cmatrix(i,9,1),i=1,3)

         call getline(fp,line,status)
         if (status.eq.1) goto 20
         read(line,*) (cmatrix(i,10,1),i=1,2)
      endif

      goto 50

 20   close(fp)
      print*,'Problem with CropParaFile',line
      stop

 50   return
      end




!soil--------------------------------------------------------------------
      subroutine getline(fileunit,line,status)

!c------------------------------------------------------------------------
      implicit none

      integer fileunit,status
      character*200 line

      status=0
 10   read(fileunit,'(a)',end=50)line
      if (line(1:1).eq.';') goto 10

      return
 50   status=1

      return
      end



!crop***********************************************************************
      subroutine initCropCommon(cmatrix,startww,emergedate, &
          sodat,laistart)
!
!     init crop specific variables
!--------------------------------------------------------------------------
      implicit none

      integer cspec,ccols,crows,nstag

      parameter (cspec=1)
      parameter (ccols=6)
      parameter (crows=10)
      parameter (nstag=5)

      common/crops/sowdate,cutday1,cutday2,cutday3,grassflag,senstage
      common/cropss/hs,conveff,minbiom,cutlai,startlai,partl,partst, &
          partrt,partgr,SLA,senfac,efv,laytck,maxrtdep,tb

      double precision tb,hs(nstag),conveff,partl(nstag),partst(nstag), &
          partrt(nstag),partgr(nstag),SLA,senfac,efv,laytck, &
          maxrtdep,minbiom,cutlai,startlai
      integer sowdate,emergedate,cutday1,cutday2,cutday3, &
          grassflag,sodat,senstage
      integer idx,i
      double precision cmatrix(ccols,crows,cspec),startww,laistart


      idx=1

      grassflag=int(cmatrix(1,1,idx)+0.5)
!e      sowdate=int(cmatrix(2,1,idx)+0.5)
      sowdate=sodat
!e      emergedate=int(cmatrix(3,1,idx)+0.5)
      emergedate=sowdate+10
      do i=1,nstag
         hs(i)=cmatrix(i,2,idx)
      enddo
      tb=cmatrix(1,3,idx)
      startww=cmatrix(2,3,idx)
      startLAI=cmatrix(3,3,idx)
      laistart=startLAI
      senstage=int(cmatrix(4,3,idx)+0.5)
      do i=1,nstag
         partl(i)=cmatrix(i,4,idx)
         partst(i)=cmatrix(i,5,idx)
         partrt(i)=cmatrix(i,6,idx)
         partgr(i)=cmatrix(i,7,idx)
      enddo
      SLA=cmatrix(1,8,idx)
      conveff=cmatrix(2,8,idx)
      efv=cmatrix(3,8,idx)
      maxrtdep=cmatrix(4,8,idx)
      laytck=cmatrix(5,8,idx)
      senfac=cmatrix(6,8,idx)
      if (grassflag.eq.1) then
         cutday1=int(cmatrix(1,9,idx)+0.5)
         cutday2=int(cmatrix(2,9,idx)+0.5)
         cutday3=int(cmatrix(3,9,idx)+0.5)
         minbiom=cmatrix(1,10,idx)
         cutlai=cmatrix(2,10,idx)
      else
         cutday1=0
         cutday2=0
         cutday3=0
         minbiom=0.
         cutlai=0.
      endif
!e      sodat=sowdate

      return
      end

